home *** CD-ROM | disk | FTP | other *** search
/ Dr. Windows 3 / dr win3.zip / dr win3 / WINPROGS / WINSRC20.ZIP / LORENZ.C < prev    next >
C/C++ Source or Header  |  1990-10-18  |  46KB  |  1,713 lines
  1. /*
  2.    This file contains two 3 dimensional orbit-type fractal
  3.    generators - IFS and LORENZ3D, along with code to generate
  4.    red/blue 3D images. Tim Wegner
  5. */
  6.  
  7. #include <stdio.h>
  8. #include <stdlib.h>
  9. #include <float.h>
  10. #include "fractint.h"
  11. #include "fractype.h"
  12. struct affine
  13. {
  14.    /* weird order so a,b,e and c,d,f are vectors */
  15.    double a;
  16.    double b;
  17.    double e;
  18.    double c;
  19.    double d;
  20.    double f;
  21. };
  22.  
  23. /* Routines in this module    */
  24.  
  25. int  orbit3dlongsetup();
  26. int  orbit3dfloatsetup();
  27. int  lorenz3dlongorbit(long *, long *, long *);
  28. int  lorenz3dfloatorbit(double *, double *, double *);
  29. int  henonfloatorbit(double *, double *, double *);
  30. int  henonlongorbit(long *, long *, long *);
  31. int  rosslerfloatorbit(double *, double *, double *);
  32. int  pickoverfloatorbit(double *, double *, double *);
  33. int  gingerbreadfloatorbit(double *, double *, double *);
  34. int  rosslerlongorbit(long *, long *, long *);
  35. int  kamtorusfloatorbit(double *, double *, double *);
  36. int  kamtoruslongorbit(long *, long *, long *);
  37. int  orbit2dfloat(void);
  38. int  orbit2dlong(void);
  39. int  orbit3dlongcalc(void);
  40. int  orbit3dfloatcalc(void);
  41. int  funny_glasses_call(int (*calc)());
  42. int  ifs(void);
  43. int  orbit3dfloat(void);
  44. int  orbit3dlong(void);
  45. int  ifs3d(void);
  46.  
  47. static int  ifs3dlong(void);
  48. static int  ifs3dfloat(void);
  49. static double determinant(double mat[3][3]);
  50. static int  solve3x3(double mat[3][3],double vec[3],double ans[3]);
  51. static int  setup_convert_to_screen(struct affine *);
  52. static void setupmatrix(MATRIX);
  53.  
  54. int realtime;
  55. extern int overflow;
  56. extern int soundflag;
  57. extern int basehertz;
  58. extern int fractype;
  59. extern int glassestype;
  60. extern int whichimage;
  61. extern int init3d[];
  62. extern char floatflag;
  63. extern VECTOR view;
  64. extern int xxadjust, yyadjust;
  65. extern int xxadjust1, yyadjust1;
  66. extern int xshift,yshift;
  67. extern int xshift1,yshift1;
  68. extern void (*plot)();
  69. extern void (*standardplot)();
  70. extern int    debugflag;            /* for debugging purposes */
  71. extern int    xdots, ydots;        /* coordinates of dots on the screen  */
  72. extern int    maxit;                /* try this many iterations */
  73. extern double param[];
  74. extern double    xxmin,xxmax,yymin,yymax,xx3rd,yy3rd; /* selected screen corners  */
  75. extern    int    diskvideo;            /* for disk-video klooges */
  76. extern int    bitshift;            /* bit shift for fudge */
  77. extern long    fudge;                /* fudge factor (2**n) */
  78. extern int    colors;             /* maximum colors available */
  79.  
  80. extern int display3d;
  81. static int t;
  82. static long l_dx,l_dy,l_dz,l_dt,l_a,l_b,l_c,l_d;
  83. static long l_adt,l_bdt,l_cdt,l_xdt,l_ydt;
  84. static long l_initx,l_inity,l_initz;
  85. static long initorbitlong[3];
  86.  
  87. static double dx,dy,dz,dt,a,b,c,d;
  88. static double adt,bdt,cdt,xdt,ydt;
  89. static double initx,inity,initz;
  90. static double initorbit[3];
  91. extern int inside;
  92.  
  93. extern int  alloc_resume(int,int);
  94. extern int  start_resume();
  95. extern void end_resume();
  96. extern int  put_resume(int len, ...);
  97. extern int  get_resume(int len, ...);
  98. extern int  calc_status, resuming;
  99.  
  100. /* these are potential user parameters */
  101. int connect = 1;    /* flag to connect points with a line */
  102. int waste = 100;    /* waste this many points before plotting */
  103. int projection = 2; /* projection plane - default is to plot x-y */
  104.  
  105. extern int active_system;
  106.  
  107. /******************************************************************/
  108. /*           zoom box conversion functions          */
  109. /******************************************************************/
  110.  
  111. static double determinant(mat) /* determinant of 3x3 matrix */
  112. double mat[3][3];
  113. {
  114.    /* calculate determinant of 3x3 matrix */
  115.    return(mat[0][0]*mat[1][1]*mat[2][2] +
  116.       mat[0][2]*mat[1][0]*mat[2][1] +
  117.       mat[0][1]*mat[1][2]*mat[2][0] -
  118.       mat[2][0]*mat[1][1]*mat[0][2] -
  119.       mat[1][0]*mat[0][1]*mat[2][2] -
  120.       mat[0][0]*mat[1][2]*mat[2][1]);
  121.  
  122. }
  123.  
  124. static int solve3x3(mat,vec,ans) /* solve 3x3 inhomogeneous linear equations */
  125. double mat[3][3], vec[3], ans[3];
  126. {
  127.    /* solve 3x3 linear equation [mat][ans] = [vec] */
  128.    double denom;
  129.    double tmp[3][3];
  130.    int i;
  131.    denom = determinant(mat);
  132.    if(fabs(denom) < DBL_EPSILON) /* test if can solve */
  133.      return(-1);
  134.    memcpy(tmp,mat,sizeof(double)*9);
  135.    for(i=0;i<3;i++)
  136.    {
  137.       tmp[0][i] = vec[0];
  138.       tmp[1][i] = vec[1];
  139.       tmp[2][i] = vec[2];
  140.       ans[i]  =  determinant(tmp)/denom;
  141.       tmp[0][i] = mat[0][i];
  142.       tmp[1][i] = mat[1][i];
  143.       tmp[2][i] = mat[2][i];
  144.     }
  145.     return(0);
  146. }
  147.  
  148.  
  149. /* Conversion of complex plane to screen coordinates for rotating zoom box.
  150.    Assume there is an affine transformation mapping complex zoom parallelogram
  151.    to rectangular screen. We know this map must map parallelogram corners to
  152.    screen corners, so we have following equations:
  153.  
  154.       a*xxmin+b*yymax+e == 0        (upper left)
  155.       c*xxmin+d*yymax+f == 0
  156.  
  157.       a*xx3rd+b*yy3rd+e == 0        (lower left)
  158.       c*xx3rd+d*yy3rd+f == ydots-1
  159.  
  160.       a*xxmax+b*yymin+e == xdots-1  (lower right)
  161.       c*xxmax+d*yymin+f == ydots-1
  162.  
  163.       First we must solve for a,b,c,d,e,f - (which we do once per image),
  164.       then we just apply the transformation to each orbit value.
  165. */
  166.  
  167. static int setup_convert_to_screen(struct affine *scrn_cnvt)
  168. {
  169.    /* we do this twice - rather than having six equations with six unknowns,
  170.       everything partitions to two sets of three equations with three
  171.       unknowns. Nice, because who wants to calculate a 6x6 determinant??
  172.     */
  173.    double mat[3][3];
  174.    double vec[3];
  175.    /*
  176.       first these equations - solve for a,b,e
  177.       a*xxmin+b*yymax+e == 0        (upper left)
  178.       a*xx3rd+b*yy3rd+e == 0        (lower left)
  179.       a*xxmax+b*yymin+e == xdots-1  (lower right)
  180.    */
  181.    mat[0][0] = xxmin;
  182.    mat[0][1] = yymax;
  183.    mat[0][2] = 1.0;
  184.    mat[1][0] = xx3rd;
  185.    mat[1][1] = yy3rd;
  186.    mat[1][2] = 1.0;
  187.    mat[2][0] = xxmax;
  188.    mat[2][1] = yymin;
  189.    mat[2][2] = 1.0;
  190.    vec[0]    = 0.0;
  191.    vec[1]    = 0.0;
  192.    vec[2]    = (double)(xdots-1);
  193.  
  194.    if(solve3x3(mat,vec, &(scrn_cnvt->a)))
  195.       return(-1);
  196.    /*
  197.       now solve these:
  198.       c*xxmin+d*yymax+f == 0
  199.       c*xx3rd+d*yy3rd+f == ydots-1
  200.       c*xxmax+d*yymin+f == ydots-1
  201.       (mat[][] has not changed - only vec[])
  202.    */
  203.    vec[0]    = 0.0;
  204.    vec[1]    = (double)(ydots-1);
  205.    vec[2]    = (double)(ydots-1);
  206.  
  207.    if(solve3x3(mat,vec, &scrn_cnvt->c))
  208.       return(-1);
  209.    return(0);
  210. }
  211.  
  212. /******************************************************************/
  213. /*   setup functions - put in fractalspecific[fractype].per_image */
  214. /******************************************************************/
  215.  
  216. double orbit;
  217. long   l_orbit;
  218.  
  219. extern double sinx,cosx;
  220. long l_sinx,l_cosx;
  221.  
  222. int orbit3dlongsetup()
  223. {
  224.    connect = 1;
  225.    waste = 100;
  226.    projection = 2;
  227.    if(fractype==LHENON || fractype==KAM || fractype==KAM3D)
  228.       connect=0;
  229.    if(fractype==LROSSLER)
  230.       waste = 500;
  231.    if(fractype==LLORENZ)
  232.       projection = 1;
  233.  
  234.    initorbitlong[0] = fudge;  /* initial conditions */
  235.    initorbitlong[1] = fudge;
  236.    initorbitlong[2] = fudge;
  237.  
  238.    if(fractype==LHENON)
  239.    {
  240.       l_a =  param[0]*fudge;
  241.       l_b =  param[1]*fudge;
  242.       l_c =  param[2]*fudge;
  243.       l_d =  param[3]*fudge;
  244.    }
  245.    else if(fractype==KAM || fractype==KAM3D)
  246.    {
  247.       a   = param[0];        /* angle */
  248.       if(param[1] <= 0.0)
  249.      param[1] = .01;
  250.       l_b =  param[1]*fudge;    /* stepsize */
  251.       l_c =  param[2]*fudge;    /* stop */
  252.       t = l_d =  param[3];     /* points per orbit */
  253.  
  254.       l_sinx = sin(a)*fudge;
  255.       l_cosx = cos(a)*fudge;
  256.       l_orbit = 0;
  257.       initorbitlong[0] = initorbitlong[1] = initorbitlong[2] = 0;
  258.    }
  259.    else
  260.    {
  261.       l_dt = param[0]*fudge;
  262.       l_a =  param[1]*fudge;
  263.       l_b =  param[2]*fudge;
  264.       l_c =  param[3]*fudge;
  265.    }
  266.  
  267.    /* precalculations for speed */
  268.    l_adt = multiply(l_a,l_dt,bitshift);
  269.    l_bdt = multiply(l_b,l_dt,bitshift);
  270.    l_cdt = multiply(l_c,l_dt,bitshift);
  271.    return(1);
  272. }
  273.  
  274.  
  275.  
  276. int orbit3dfloatsetup()
  277. {
  278.    connect = 1;
  279.    waste = 100;
  280.    projection = 2;
  281.  
  282.    if(fractype==FPHENON || fractype==FPPICKOVER || fractype==FPGINGERBREAD
  283.         || fractype == KAMFP || fractype == KA